load data, libraries, functions

library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(MOFA2)
## Warning: package 'MOFA2' was built under R version 4.3.2
## 
## Attaching package: 'MOFA2'
## 
## The following object is masked from 'package:stats':
## 
##     predict
library(ggplot2)
library(MOFAdata)

Day0 T cell pol

Age

Factor 15

MOFAobject = load_model('../data/MOFA_models/MOFA2_noScale_train2020_21.hdf5')
MOFAobject
## Trained MOFA with the following characteristics: 
##  Number of views: 6 
##  Views names: geneExp ab cytokineL cell_freq tcellPol tcellAct 
##  Number of features (per view): 6650 31 14 39 1 2 
##  Number of groups: 1 
##  Groups names: group1 
##  Number of samples (per group): 158 
##  Number of factors: 15
plot_variance_explained(MOFAobject, x="view", y="factor")

plot_factor(MOFAobject, 
            factor = 15,
            color_by = "IgG_PT"
)

plot_weights(MOFAobject,
             #view = "geneExp",
             factor = 15,
             nfeatures = 20,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "geneExp",
             factor = 15,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "cell_freq",
             factor = 15,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "ab",
             factor = 15,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "cytokineL",
             factor = 15,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_top_weights(MOFAobject,
                 #view = "geneExp",
                 factor = 15,
                 nfeatures = 40
)

plot_top_weights(MOFAobject,
                 view = "cell_freq",
                 factor = 15,
                 nfeatures = 40
)

plot_data_heatmap(MOFAobject,
                  view = "geneExp",         # view of interest
                  factor = 15,             # factor of interest
                  features = 60,          # number of features to plot (they are selected by weight)
                  
                  # extra arguments that are passed to the `pheatmap` function
                  cluster_rows = TRUE, cluster_cols = FALSE,
                  show_rownames = TRUE, show_colnames = FALSE
)

plot_data_scatter(MOFAobject,
                  view = "geneExp",         # view of interest
                  factor = 15,             # factor of interest
                  features = 11,           # number of features to plot (they are selected by weight)
                  add_lm = TRUE,         # add linear regression
                  color_by = "Meta.infancy_vac"
)

plot_data_scatter(MOFAobject,
                  view = "cytokineL",         # view of interest
                  factor = 15,             # factor of interest
                  features = 11,           # number of features to plot (they are selected by weight)
                  add_lm = TRUE,         # add linear regression
                  color_by = "Meta.infancy_vac"
)

factors <- get_factors(MOFAobject, factors = "all")
lapply(factors,dim)
## $group1
## [1] 158  15
weights <- get_weights(MOFAobject, views = "all", factors = "Factor15")
a= weights$geneExp
ggplot(a, aes(x=Factor15)) +
    geom_density() 

a['Expr.ENSG00000277632.1', ]
## [1] -0.0363342
a = a[order(a[,1]), ]
head(a)
## Expr.ENSG00000166508.17 Expr.ENSG00000038219.12 Expr.ENSG00000141027.20 
##              -0.3850143              -0.2863488              -0.2863296 
##  Expr.ENSG00000228474.5 Expr.ENSG00000070814.17 Expr.ENSG00000166046.10 
##              -0.2683452              -0.2596098              -0.2537850
tail(a)
## Expr.ENSG00000135002.11  Expr.ENSG00000164713.9 Expr.ENSG00000059378.12 
##               0.4043662               0.4137762               0.4394598 
## Expr.ENSG00000134970.13 Expr.ENSG00000005812.10  Expr.ENSG00000126821.7 
##               0.4603818               0.4899954               0.5704848
tail(names(a), n=200)
##   [1] "Expr.ENSG00000111269.2"  "Expr.ENSG00000106785.14"
##   [3] "Expr.ENSG00000164164.15" "Expr.ENSG00000105963.13"
##   [5] "Expr.ENSG00000148572.15" "Expr.ENSG00000095906.16"
##   [7] "Expr.ENSG00000198728.10" "Expr.ENSG00000118961.14"
##   [9] "Expr.ENSG00000161091.12" "Expr.ENSG00000170836.11"
##  [11] "Expr.ENSG00000176393.10" "Expr.ENSG00000137841.11"
##  [13] "Expr.ENSG00000185127.6"  "Expr.ENSG00000185947.14"
##  [15] "Expr.ENSG00000244462.7"  "Expr.ENSG00000147439.12"
##  [17] "Expr.ENSG00000117090.14" "Expr.ENSG00000143337.18"
##  [19] "Expr.ENSG00000140511.11" "Expr.ENSG00000213064.9" 
##  [21] "Expr.ENSG00000118596.11" "Expr.ENSG00000178951.8" 
##  [23] "Expr.ENSG00000166822.12" "Expr.ENSG00000146232.15"
##  [25] "Expr.ENSG00000100106.19" "Expr.ENSG00000120925.15"
##  [27] "Expr.ENSG00000102100.15" "Expr.ENSG00000167747.14"
##  [29] "Expr.ENSG00000008282.8"  "Expr.ENSG00000110172.11"
##  [31] "Expr.ENSG00000185728.16" "Expr.ENSG00000164933.11"
##  [33] "Expr.ENSG00000130560.8"  "Expr.ENSG00000051523.10"
##  [35] "Expr.ENSG00000142065.13" "Expr.ENSG00000109445.10"
##  [37] "Expr.ENSG00000183513.8"  "Expr.ENSG00000134215.15"
##  [39] "Expr.ENSG00000165156.14" "Expr.ENSG00000170619.9" 
##  [41] "Expr.ENSG00000139826.5"  "Expr.ENSG00000180530.10"
##  [43] "Expr.ENSG00000163527.9"  "Expr.ENSG00000085365.17"
##  [45] "Expr.ENSG00000114450.9"  "Expr.ENSG00000132670.20"
##  [47] "Expr.ENSG00000114735.9"  "Expr.ENSG00000125826.20"
##  [49] "Expr.ENSG00000181904.8"  "Expr.ENSG00000148110.15"
##  [51] "Expr.ENSG00000153814.11" "Expr.ENSG00000053918.15"
##  [53] "Expr.ENSG00000198382.8"  "Expr.ENSG00000160446.18"
##  [55] "Expr.ENSG00000134851.12" "Expr.ENSG00000086300.15"
##  [57] "Expr.ENSG00000182325.10" "Expr.ENSG00000198176.12"
##  [59] "Expr.ENSG00000276045.2"  "Expr.ENSG00000025708.13"
##  [61] "Expr.ENSG00000167543.15" "Expr.ENSG00000154174.7" 
##  [63] "Expr.ENSG00000132286.11" "Expr.ENSG00000205268.10"
##  [65] "Expr.ENSG00000165458.13" "Expr.ENSG00000254087.7" 
##  [67] "Expr.ENSG00000116044.15" "Expr.ENSG00000170779.10"
##  [69] "Expr.ENSG00000109270.12" "Expr.ENSG00000196646.11"
##  [71] "Expr.ENSG00000106993.11" "Expr.ENSG00000106460.18"
##  [73] "Expr.ENSG00000111912.19" "Expr.ENSG00000135597.18"
##  [75] "Expr.ENSG00000221869.4"  "Expr.ENSG00000133703.11"
##  [77] "Expr.ENSG00000071994.10" "Expr.ENSG00000157837.15"
##  [79] "Expr.ENSG00000111875.7"  "Expr.ENSG00000171612.6" 
##  [81] "Expr.ENSG00000188452.13" "Expr.ENSG00000059758.7" 
##  [83] "Expr.ENSG00000099849.14" "Expr.ENSG00000075568.16"
##  [85] "Expr.ENSG00000168159.11" "Expr.ENSG00000084093.16"
##  [87] "Expr.ENSG00000125520.13" "Expr.ENSG00000106392.10"
##  [89] "Expr.ENSG00000104518.10" "Expr.ENSG00000184990.12"
##  [91] "Expr.ENSG00000070476.14" "Expr.ENSG00000165476.13"
##  [93] "Expr.ENSG00000167208.14" "Expr.ENSG00000163683.11"
##  [95] "Expr.ENSG00000133997.11" "Expr.ENSG00000151135.9" 
##  [97] "Expr.ENSG00000161904.11" "Expr.ENSG00000140848.16"
##  [99] "Expr.ENSG00000117906.13" "Expr.ENSG00000121774.17"
## [101] "Expr.ENSG00000175416.12" "Expr.ENSG00000103266.10"
## [103] "Expr.ENSG00000178562.17" "Expr.ENSG00000185000.11"
## [105] "Expr.ENSG00000103187.7"  "Expr.ENSG00000163344.5" 
## [107] "Expr.ENSG00000105135.15" "Expr.ENSG00000106263.17"
## [109] "Expr.ENSG00000185896.10" "Expr.ENSG00000130734.9" 
## [111] "Expr.ENSG00000188811.13" "Expr.ENSG00000165782.10"
## [113] "Expr.ENSG00000104852.14" "Expr.ENSG00000185024.16"
## [115] "Expr.ENSG00000048405.9"  "Expr.ENSG00000162783.10"
## [117] "Expr.ENSG00000100242.15" "Expr.ENSG00000069399.14"
## [119] "Expr.ENSG00000189319.13" "Expr.ENSG00000142751.14"
## [121] "Expr.ENSG00000105447.12" "Expr.ENSG00000164305.18"
## [123] "Expr.ENSG00000181788.3"  "Expr.ENSG00000138231.12"
## [125] "Expr.ENSG00000161677.11" "Expr.ENSG00000171488.14"
## [127] "Expr.ENSG00000159873.9"  "Expr.ENSG00000143514.16"
## [129] "Expr.ENSG00000196449.3"  "Expr.ENSG00000254470.2" 
## [131] "Expr.ENSG00000138032.20" "Expr.ENSG00000172164.13"
## [133] "Expr.ENSG00000159363.17" "Expr.ENSG00000184897.5" 
## [135] "Expr.ENSG00000130303.12" "Expr.ENSG00000108604.15"
## [137] "Expr.ENSG00000178605.13" "Expr.ENSG00000130511.15"
## [139] "Expr.ENSG00000196843.15" "Expr.ENSG00000172216.5" 
## [141] "Expr.ENSG00000000419.12" "Expr.ENSG00000116649.9" 
## [143] "Expr.ENSG00000198455.4"  "Expr.ENSG00000140367.11"
## [145] "Expr.ENSG00000119900.7"  "Expr.ENSG00000171603.16"
## [147] "Expr.ENSG00000072401.14" "Expr.ENSG00000198087.7" 
## [149] "Expr.ENSG00000108784.9"  "Expr.ENSG00000198055.10"
## [151] "Expr.ENSG00000122490.18" "Expr.ENSG00000134758.13"
## [153] "Expr.ENSG00000104671.7"  "Expr.ENSG00000165802.22"
## [155] "Expr.ENSG00000139977.13" "Expr.ENSG00000198133.8" 
## [157] "Expr.ENSG00000168056.15" "Expr.ENSG00000140280.13"
## [159] "Expr.ENSG00000065000.15" "Expr.ENSG00000099624.7" 
## [161] "Expr.ENSG00000166128.12" "Expr.ENSG00000173598.13"
## [163] "Expr.ENSG00000108679.12" "Expr.ENSG00000184207.8" 
## [165] "Expr.ENSG00000148341.17" "Expr.ENSG00000168071.21"
## [167] "Expr.ENSG00000174738.12" "Expr.ENSG00000149658.17"
## [169] "Expr.ENSG00000071859.14" "Expr.ENSG00000146425.10"
## [171] "Expr.ENSG00000185129.5"  "Expr.ENSG00000135363.11"
## [173] "Expr.ENSG00000111711.9"  "Expr.ENSG00000066044.14"
## [175] "Expr.ENSG00000111666.10" "Expr.ENSG00000080371.5" 
## [177] "Expr.ENSG00000125910.5"  "Expr.ENSG00000136238.17"
## [179] "Expr.ENSG00000169660.15" "Expr.ENSG00000105122.12"
## [181] "Expr.ENSG00000151553.14" "Expr.ENSG00000148943.11"
## [183] "Expr.ENSG00000170385.9"  "Expr.ENSG00000099804.8" 
## [185] "Expr.ENSG00000175467.14" "Expr.ENSG00000166337.9" 
## [187] "Expr.ENSG00000136897.7"  "Expr.ENSG00000131979.18"
## [189] "Expr.ENSG00000060491.16" "Expr.ENSG00000162368.13"
## [191] "Expr.ENSG00000259956.1"  "Expr.ENSG00000116863.10"
## [193] "Expr.ENSG00000135272.10" "Expr.ENSG00000123933.16"
## [195] "Expr.ENSG00000135002.11" "Expr.ENSG00000164713.9" 
## [197] "Expr.ENSG00000059378.12" "Expr.ENSG00000134970.13"
## [199] "Expr.ENSG00000005812.10" "Expr.ENSG00000126821.7"
b= weights$ab
ggplot(b, aes(x=Factor15)) +
    geom_density() 

c= weights$cytokineL
ggplot(c, aes(x=Factor15)) +
    geom_density() 

d= weights$cell_freq
ggplot(d, aes(x=Factor15)) +
    geom_density() 

enrichment

against Hallmark

MOFAobject3 = MOFAobject
ens_map = read_tsv('~/Documents/Projects/CMI-PB2/CMI-PB-Oct2023-FINAL/data/raw_prediction_dataset/gene_90_38_export.tsv')
## Rows: 63971 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): biotype, versioned_ensembl_gene_id, description, display_label, syn...
## dbl (6): seq_region_strand, seq_region_start, seq_region_end, canonical_tran...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(ens_map)
## # A tibble: 6 × 12
##   biotype        seq_region_strand seq_region_start seq_region_end
##   <chr>                      <dbl>            <dbl>          <dbl>
## 1 miRNA                         -1         55372940       55373034
## 2 snRNA                         -1         59450673       59450772
## 3 miRNA                          1         41691585       41691678
## 4 protein_coding                 1         54201473       54208260
## 5 misc_RNA                       1         10287413       10287482
## 6 snRNA                          1        109992325      109992431
## # ℹ 8 more variables: versioned_ensembl_gene_id <chr>, description <chr>,
## #   canonical_transcript_id <dbl>, display_label <chr>, synonym <chr>,
## #   dbprimary_acc <dbl>, gene_id <dbl>, gene_summary <chr>
genes = as.character(rownames(MOFAobject3@data$geneExp$group1))
genes = gsub('Expr.', '', genes)
table(genes %in% ens_map$versioned_ensembl_gene_id)
## 
## TRUE 
## 6650
ens_map = ens_map[!duplicated(ens_map$versioned_ensembl_gene_id), ]
rownames(ens_map) = ens_map$versioned_ensembl_gene_id
## Warning: Setting row names on a tibble is deprecated.
genes_sym = ens_map[genes, 'display_label']$display_label

rownames(MOFAobject3@data$geneExp$group1) = genes_sym
names(MOFAobject3@intercepts$geneExp$group1) = genes_sym

MOFAobject3@features_metadata$feature[1:6650] = genes_sym

rownames(MOFAobject3@expectations$W$geneExp) = genes_sym
hallmark = fgsea::gmtPathways('~/Documents/References/Genesets/MSigDB/h.all.v2023.2.Hs.symbols.gmt.txt')


# Get all unique genes
all_genes <- sort(unique(unlist(hallmark)))

# Create the binary matrix
binary_matrix <- sapply(all_genes, function(g) {
    as.numeric(sapply(hallmark, function(x) g %in% x))
})

# Transpose and convert to a data frame for better visualization
#binary_matrix <- as.data.frame(t(binary_matrix))
binary_matrix[1:5, 1:5]
##      A2M AAAS AADAT AARS1 ABAT
## [1,]   0    0     0     0    0
## [2,]   0    0     0     1    0
## [3,]   0    0     0     0    0
## [4,]   0    0     0     0    0
## [5,]   0    0     0     0    0
# Add row names
rownames(binary_matrix) <- names(hallmark)

# Display the binary matrix
binary_matrix[1:5, 1:10]
##                              A2M AAAS AADAT AARS1 ABAT ABCA1 ABCA2 ABCA3 ABCA4
## HALLMARK_ADIPOGENESIS          0    0     0     0    0     1     0     0     0
## HALLMARK_ALLOGRAFT_REJECTION   0    0     0     1    0     0     0     0     0
## HALLMARK_ANDROGEN_RESPONSE     0    0     0     0    0     0     0     0     0
## HALLMARK_ANGIOGENESIS          0    0     0     0    0     0     0     0     0
## HALLMARK_APICAL_JUNCTION       0    0     0     0    0     0     0     0     0
##                              ABCA5
## HALLMARK_ADIPOGENESIS            0
## HALLMARK_ALLOGRAFT_REJECTION     0
## HALLMARK_ANDROGEN_RESPONSE       0
## HALLMARK_ANGIOGENESIS            0
## HALLMARK_APICAL_JUNCTION         0
dim(binary_matrix)
## [1]   50 4384
enrichment3.parametric.pos <- run_enrichment(MOFAobject3,
                                            view = "geneExp", factors = 15,
                                            feature.sets = as.matrix(binary_matrix),
                                            sign = "positive",
                                            statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 2153 features.
## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 48 
## Set statistic: mean.diff 
## Statistical test: parametric
## Subsetting weights with positive sign
## 
plot_enrichment(enrichment3.parametric.pos, 
                factor = 'Factor15', 
                max.pathways = 15
)

plot_enrichment_detailed(enrichment3.parametric.pos, 
                         factor = 'Factor15', 
                         max.genes = 8, 
                         max.pathways = 5
)

enrichment3.parametric.neg <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = 15,
                                             feature.sets = as.matrix(binary_matrix),
                                             sign = "negative",
                                             statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 2153 features.
## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 48 
## Set statistic: mean.diff 
## Statistical test: parametric
## Subsetting weights with negative sign
## 

against BTM

btm = fgsea::gmtPathways('~/Documents/References/Genesets/BTM/BTM_for_GSEA_20131008.gmt')
lapply(btm, length)
## $`targets of FOSL1/2 (M0)`
## [1] 12
## 
## $`integrin cell surface interactions (I) (M1.0)`
## [1] 29
## 
## $`integrin cell surface interactions (II) (M1.1)`
## [1] 12
## 
## $`extracellular matrix (I) (M2.0)`
## [1] 30
## 
## $`extracellular matrix (II) (M2.1)`
## [1] 45
## 
## $`extracellular matrix (III) (M2.2)`
## [1] 14
## 
## $`regulation of signal transduction (M3)`
## [1] 47
## 
## $`cell cycle and transcription (M4.0)`
## [1] 335
## 
## $`cell cycle (I) (M4.1)`
## [1] 145
## 
## $`PLK1 signaling events (M4.2)`
## [1] 34
## 
## $`myeloid cell enriched receptors and transporters (M4.3)`
## [1] 31
## 
## $`mitotic cell cycle - DNA replication (M4.4)`
## [1] 30
## 
## $`mitotic cell cycle in stimulated CD4 T cells (M4.5)`
## [1] 35
## 
## $`cell division in stimulated CD4 T cells (M4.6)`
## [1] 20
## 
## $`mitotic cell cycle (M4.7)`
## [1] 21
## 
## $`cell division - E2F transcription network (M4.8)`
## [1] 19
## 
## $`mitotic cell cycle in stimulated CD4 T cells (M4.9)`
## [1] 16
## 
## $`cell cycle (II) (M4.10)`
## [1] 14
## 
## $`mitotic cell cycle in stimulated CD4 T cells (M4.11)`
## [1] 12
## 
## $`C-MYC transcriptional network (M4.12)`
## [1] 12
## 
## $`cell junction (GO) (M4.13)`
## [1] 11
## 
## $`Rho GTPase cycle (M4.14)`
## [1] 10
## 
## $`enriched in monocytes (I) (M4.15)`
## [1] 11
## 
## $`regulation of antigen presentation and immune response (M5.0)`
## [1] 81
## 
## $`T cell activation and signaling (M5.1)`
## [1] 25
## 
## $`mitotic cell division (M6)`
## [1] 32
## 
## $`enriched in T cells (I) (M7.0)`
## [1] 62
## 
## $`T cell activation (I) (M7.1)`
## [1] 50
## 
## $`enriched in NK cells (I) (M7.2)`
## [1] 49
## 
## $`T cell activation (II) (M7.3)`
## [1] 31
## 
## $`T cell activation (III) (M7.4)`
## [1] 14
## 
## $`E2F transcription factor network (M8)`
## [1] 18
## 
## $`B cell development (M9)`
## [1] 11
## 
## $`E2F1 targets (Q3) (M10.0)`
## [1] 33
## 
## $`E2F1 targets (Q4) (M10.1)`
## [1] 21
## 
## $`enriched in monocytes (II) (M11.0)`
## [1] 189
## 
## $`blood coagulation (M11.1)`
## [1] 22
## 
## $`formyl peptide receptor mediated neutrophil response (M11.2)`
## [1] 10
## 
## $`CD28 costimulation (M12)`
## [1] 10
## 
## $`innate activation by cytosolic DNA sensing (M13)`
## [1] 11
## 
## $`T cell differentiation (M14)`
## [1] 12
## 
## $`Ran mediated mitosis (M15)`
## [1] 13
## 
## $`TLR and inflammatory signaling (M16)`
## [1] 45
## 
## $`Hox cluster I (M17.0)`
## [1] 16
## 
## $`Hox cluster II (M17.1)`
## [1] 12
## 
## $`Hox cluster III (M17.2)`
## [1] 12
## 
## $`Hox cluster IV (M17.3)`
## [1] 10
## 
## $`T cell differentiation via ITK and PKC (M18)`
## [1] 11
## 
## $`T cell differentiation (Th2) (M19)`
## [1] 17
## 
## $`AP-1 transcription factor network (M20)`
## [1] 15
## 
## $`cell adhesion (lymphocyte homing) (M21)`
## [1] 10
## 
## $`mismatch repair (I) (M22.0)`
## [1] 27
## 
## $`mismatch repair (II) (M22.1)`
## [1] 13
## 
## $`RA, WNT, CSF receptors network (monocyte) (M23)`
## [1] 12
## 
## $`cell activation (IL15, IL23, TNF) (M24)`
## [1] 15
## 
## $`TLR8-BAFF network (M25)`
## [1] 10
## 
## $`TBA (M26.0)`
## [1] 30
## 
## $`TBA (M26.1)`
## [1] 22
## 
## $`TBA (M26.2)`
## [1] 14
## 
## $`chemokine cluster (I) (M27.0)`
## [1] 26
## 
## $`chemokine cluster (II) (M27.1)`
## [1] 15
## 
## $`antigen presentation (lipids and proteins) (M28)`
## [1] 11
## 
## $`proinflammatory cytokines and chemokines (M29)`
## [1] 10
## 
## $`cell movement, Adhesion & Platelet activation (M30)`
## [1] 20
## 
## $`cell cycle and growth arrest (M31)`
## [1] 12
## 
## $`platelet activation (I) (M32.0)`
## [1] 23
## 
## $`platelet activation (II) (M32.1)`
## [1] 21
## 
## $`CORO1A-DEF6 network (I) (M32.2)`
## [1] 20
## 
## $`KLF12 targets network (M32.3)`
## [1] 17
## 
## $`CORO1A-DEF6 network (II) (M32.4)`
## [1] 15
## 
## $`TBA (M32.5)`
## [1] 15
## 
## $`TBA (M32.6)`
## [1] 13
## 
## $`TBA (M32.7)`
## [1] 11
## 
## $`cytoskeletal remodeling (M32.8)`
## [1] 10
## 
## $`inflammatory response (M33)`
## [1] 11
## 
## $`cytoskeletal remodeling (enriched for SRF targets) (M34)`
## [1] 10
## 
## $`signaling in T cells (I) (M35.0)`
## [1] 14
## 
## $`signaling in T cells (II) (M35.1)`
## [1] 11
## 
## $`T cell surface, activation (M36)`
## [1] 12
## 
## $`immune activation - generic cluster (M37.0)`
## [1] 347
## 
## $`enriched in neutrophils (I) (M37.1)`
## [1] 52
## 
## $`endoplasmic reticulum (M37.2)`
## [1] 19
## 
## $`cell division (M37.3)`
## [1] 10
## 
## $`chemokines and receptors (M38)`
## [1] 11
## 
## $`integrin mediated leukocyte migration (M39)`
## [1] 10
## 
## $`complement and other receptors in DCs (M40)`
## [1] 13
## 
## $`TBA (M41.0)`
## [1] 17
## 
## $`TBA (M41.1)`
## [1] 15
## 
## $`TBA (M41.2)`
## [1] 14
## 
## $`TBA (M41.3)`
## [1] 10
## 
## $`ATF targets network (M41.4)`
## [1] 10
## 
## $`platelet activation (III) (M42)`
## [1] 10
## 
## $`myeloid, dendritic cell activation via NFkB (I) (M43.0)`
## [1] 15
## 
## $`myeloid, dendritic cell activation via NFkB (II) (M43.1)`
## [1] 14
## 
## $`T cell signaling and costimulation (M44)`
## [1] 10
## 
## $`leukocyte activation and migration (M45)`
## [1] 11
## 
## $`cell division (stimulated CD4+ T cells) (M46)`
## [1] 28
## 
## $`enriched in B cells (I) (M47.0)`
## [1] 53
## 
## $`enriched in B cells (II) (M47.1)`
## [1] 43
## 
## $`enriched in B cells (III) (M47.2)`
## [1] 37
## 
## $`enriched in B cells (IV) (M47.3)`
## [1] 16
## 
## $`enriched in B cells (V) (M47.4)`
## [1] 11
## 
## $`TBA (M48)`
## [1] 13
## 
## $`transcription regulation in cell development (M49)`
## [1] 47
## 
## $`CD1 and other DC receptors (M50)`
## [1] 10
## 
## $`cell adhesion (M51)`
## [1] 38
## 
## $`T cell activation (IV) (M52)`
## [1] 13
## 
## $`inflammasome receptors and signaling (M53)`
## [1] 12
## 
## $`BCR signaling (M54)`
## [1] 12
## 
## $`TBA (M55)`
## [1] 12
## 
## $`suppression of MAPK signaling (M56)`
## [1] 12
## 
## $`immuregulation - monocytes, T and B cells (M57)`
## [1] 13
## 
## $`B cell development/activation (M58)`
## [1] 11
## 
## $`CCR1, 7 and cell signaling (M59)`
## [1] 10
## 
## $`lymphocyte generic cluster (M60)`
## [1] 17
## 
## $`enriched in NK cells (II) (M61.0)`
## [1] 24
## 
## $`enriched in NK cells (KIR cluster) (M61.1)`
## [1] 13
## 
## $`enriched in NK cells (receptor activation) (M61.2)`
## [1] 16
## 
## $`T & B cell development, activation (M62.0)`
## [1] 44
## 
## $`enriched for unknown TF motif CTCNANGTGNY (M62.1)`
## [1] 12
## 
## $`regulation of localization (GO) (M63)`
## [1] 12
## 
## $`enriched in activated dendritic cells/monocytes (M64)`
## [1] 17
## 
## $`IL2, IL7, TCR network (M65)`
## [1] 12
## 
## $`TBA (M66)`
## [1] 17
## 
## $`activated dendritic cells (M67)`
## [1] 12
## 
## $`RIG-1 like receptor signaling (M68)`
## [1] 10
## 
## $`enriched in B cells (VI) (M69)`
## [1] 20
## 
## $`TBA (M70.0)`
## [1] 20
## 
## $`TBA (M70.1)`
## [1] 10
## 
## $`enriched in antigen presentation (I) (M71)`
## [1] 18
## 
## $`TBA (M72.0)`
## [1] 24
## 
## $`TBA (M72.1)`
## [1] 17
## 
## $`TBA (M72.2)`
## [1] 12
## 
## $`enriched in monocytes (III) (M73)`
## [1] 12
## 
## $`transcriptional targets of glucocorticoid receptor (M74)`
## [1] 12
## 
## $`antiviral IFN signature (M75)`
## [1] 22
## 
## $`DNA repair (M76)`
## [1] 22
## 
## $`collagen, TGFB family et al (M77)`
## [1] 35
## 
## $`myeloid cell cytokines, metallopeptidases and laminins (M78)`
## [1] 11
## 
## $`TBA (M79)`
## [1] 10
## 
## $`TBA (M80)`
## [1] 20
## 
## $`enriched in myeloid cells and monocytes (M81)`
## [1] 35
## 
## $`signal transduction, plasma membrane (M82)`
## [1] 13
## 
## $`enriched in naive and memory B cells (M83)`
## [1] 10
## 
## $`integrins and cell adhesion (M84)`
## [1] 10
## 
## $`platelet activation and degranulation (M85)`
## [1] 12
## 
## $`chemokines and inflammatory molecules in myeloid cells (M86.0)`
## [1] 19
## 
## $`proinflammatory dendritic cell, myeloid cell response (M86.1)`
## [1] 13
## 
## $`transmembrane transport (I) (M87)`
## [1] 24
## 
## $`leukocyte migration (M88.0)`
## [1] 51
## 
## $`enriched in hepatocyte nuclear factors (I) (M88.1)`
## [1] 13
## 
## $`enriched in hepatocyte nuclear factors (II) (M88.2)`
## [1] 12
## 
## $`putative targets of PAX3 (M89.0)`
## [1] 16
## 
## $`putative targets of PAX3 (M89.1)`
## [1] 10
## 
## $`TBA (M90)`
## [1] 12
## 
## $`adhesion and migration, chemotaxis (M91)`
## [1] 13
## 
## $`lipid metabolism, endoplasmic reticulum (M92)`
## [1] 10
## 
## $`TBA (M93)`
## [1] 10
## 
## $`growth factor induced, enriched in nuclear receptor subfamily 4 (M94)`
## [1] 12
## 
## $`enriched in antigen presentation (II) (M95.0)`
## [1] 24
## 
## $`enriched in antigen presentation (III) (M95.1)`
## [1] 15
## 
## $`Hox cluster V (M96)`
## [1] 10
## 
## $`enriched for SMAD2/3 signaling (M97)`
## [1] 10
## 
## $`TBA (M98.0)`
## [1] 18
## 
## $`TBA (M98.1)`
## [1] 10
## 
## $`TBA (M99)`
## [1] 10
## 
## $`MAPK, RAS signaling (M100)`
## [1] 10
## 
## $`phosphatidylinositol signaling system (M101)`
## [1] 14
## 
## $`TBA (M102)`
## [1] 12
## 
## $`cell cycle (III) (M103)`
## [1] 51
## 
## $`TBA (M104)`
## [1] 10
## 
## $`TBA (M105)`
## [1] 18
## 
## $`nuclear pore complex (M106.0)`
## [1] 17
## 
## $`nuclear pore complex (mitosis) (M106.1)`
## [1] 11
## 
## $`Hox cluster VI (M107)`
## [1] 11
## 
## $`TBA (M108)`
## [1] 11
## 
## $`receptors, cell migration (M109)`
## [1] 15
## 
## $`axon guidance (M110)`
## [1] 10
## 
## $`viral sensing & immunity; IRF2 targets network (I) (M111.0)`
## [1] 17
## 
## $`viral sensing & immunity; IRF2 targets network (II) (M111.1)`
## [1] 11
## 
## $`complement activation (I) (M112.0)`
## [1] 17
## 
## $`complement activation (II) (M112.1)`
## [1] 10
## 
## $`golgi membrane (I) (M113)`
## [1] 10
## 
## $`TBA (M114.0)`
## [1] 39
## 
## $`glycerophospholipid metabolism (M114.1)`
## [1] 11
## 
## $`cytokines - recepters cluster (M115)`
## [1] 11
## 
## $`TBA (M116)`
## [1] 13
## 
## $`cell adhesion (GO) (M117)`
## [1] 21
## 
## $`enriched in monocytes (IV) (M118.0)`
## [1] 57
## 
## $`enriched in monocytes (surface) (M118.1)`
## [1] 17
## 
## $`enriched in activated dendritic cells (I) (M119)`
## [1] 11
## 
## $`TBA (M120)`
## [1] 11
## 
## $`TBA (M121)`
## [1] 12
## 
## $`enriched for cell migration (M122)`
## [1] 11
## 
## $`enriched in B cell differentiation (M123)`
## [1] 13
## 
## $`enriched in membrane proteins (M124)`
## [1] 19
## 
## $`TBA (M125)`
## [1] 11
## 
## $`double positive thymocytes (M126)`
## [1] 13
## 
## $`type I interferon response (M127)`
## [1] 12
## 
## $`TBA (M128)`
## [1] 12
## 
## $`inositol phosphate metabolism (M129)`
## [1] 13
## 
## $`enriched in G-protein coupled receptors (M130)`
## [1] 10
## 
## $`TBA (M131)`
## [1] 14
## 
## $`recruitment of neutrophils (M132)`
## [1] 11
## 
## $`cell adhesion, membrane (M133.0)`
## [1] 14
## 
## $`cell cell adhesion (M133.1)`
## [1] 11
## 
## $`Membrane, ER proteins (M134)`
## [1] 17
## 
## $`enriched in plasma membrane proteins (I) (M135.0)`
## [1] 16
## 
## $`enriched in plasma membrane proteins (II) (M135.1)`
## [1] 15
## 
## $`TBA (M136)`
## [1] 17
## 
## $`TBA (M137)`
## [1] 16
## 
## $`enriched for ubiquitination (M138)`
## [1] 11
## 
## $`lysosomal/endosomal proteins (M139)`
## [1] 11
## 
## $`extracellular matrix, complement (M140)`
## [1] 14
## 
## $`TBA (M141)`
## [1] 27
## 
## $`transmembrane and ion transporters (I) (M142)`
## [1] 13
## 
## $`nuclear pore, transport; mRNA splicing, processing (M143)`
## [1] 11
## 
## $`cell cycle, ATP binding (M144)`
## [1] 16
## 
## $`cytoskeleton/actin (SRF transcription targets) (M145.0)`
## [1] 16
## 
## $`cytoskeleton/actin (SRF transcription targets) (M145.1)`
## [1] 14
## 
## $`MHC-TLR7-TLR8 cluster (M146)`
## [1] 17
## 
## $`intracellular transport (M147)`
## [1] 17
## 
## $`TBA (M148)`
## [1] 11
## 
## $`TBA (M149)`
## [1] 10
## 
## $`innate antiviral response (M150)`
## [1] 12
## 
## $`TBA (M151)`
## [1] 10
## 
## $`TBA (source: B cells) (M152.0)`
## [1] 25
## 
## $`TBA (source: naive B cells) (M152.1)`
## [1] 21
## 
## $`TBA (source: memory B cells) (M152.2)`
## [1] 18
## 
## $`TBA (M153)`
## [1] 16
## 
## $`amino acid metabolishm and transport (M154.0)`
## [1] 33
## 
## $`transmembrane transport (SLC cluster) (M154.1)`
## [1] 17
## 
## $`G protein coupled receptors cluster (M155)`
## [1] 10
## 
## $`plasma cells & B cells, immunoglobulins (M156.0)`
## [1] 43
## 
## $`plasma cells, immunoglobulins (M156.1)`
## [1] 36
## 
## $`enriched in NK cells (III) (M157)`
## [1] 14
## 
## $`interferon alpha response (I) (M158.0)`
## [1] 16
## 
## $`interferon alpha response (II) (M158.1)`
## [1] 13
## 
## $`G protein mediated calcium signaling (M159)`
## [1] 10
## 
## $`leukocyte differentiation (M160)`
## [1] 16
## 
## $`TBA (M161)`
## [1] 14
## 
## $`plasma membrane, cell junction (M162.0)`
## [1] 18
## 
## $`cell junction (M162.1)`
## [1] 11
## 
## $`enriched in neutrophils (II) (M163)`
## [1] 14
## 
## $`xenobiotic metabolism (M164)`
## [1] 10
## 
## $`enriched in activated dendritic cells (II) (M165)`
## [1] 35
## 
## $`TBA (M166)`
## [1] 17
## 
## $`enriched in cell cycle (M167)`
## [1] 15
## 
## $`enriched in dendritic cells (M168)`
## [1] 19
## 
## $`mitosis (TF motif CCAATNNSNNNGCG) (M169)`
## [1] 16
## 
## $`TBA (M170)`
## [1] 12
## 
## $`heme biosynthesis (I) (M171)`
## [1] 11
## 
## $`enriched for TF motif TTCNRGNNNNTTC (M172)`
## [1] 10
## 
## $`erythrocyte differentiation (M173)`
## [1] 14
## 
## $`TBA (M174)`
## [1] 24
## 
## $`cell development (M175)`
## [1] 13
## 
## $`TBA (M176)`
## [1] 10
## 
## $`TBA (M177.0)`
## [1] 34
## 
## $`TBA (M177.1)`
## [1] 15
## 
## $`enriched for promoter motif NATCACGTGAY (putative SREBF1 targets) (M178)`
## [1] 10
## 
## $`enriched for TF motif PAX3 (M179)`
## [1] 10
## 
## $`TBA (M180)`
## [1] 11
## 
## $`nucleotide metabolism (M181)`
## [1] 10
## 
## $`enriched in DNA interacting proteins (M182)`
## [1] 16
## 
## $`TBA (M183)`
## [1] 10
## 
## $`TBA (M184.0)`
## [1] 18
## 
## $`TBA (M184.1)`
## [1] 11
## 
## $`TBA (M185)`
## [1] 14
## 
## $`TBA (M186)`
## [1] 11
## 
## $`metabolism in mitochondria, peroxisome (M187)`
## [1] 12
## 
## $`TBA (M188)`
## [1] 10
## 
## $`extracellular region cluster (GO) (M189)`
## [1] 16
## 
## $`TBA (M190)`
## [1] 11
## 
## $`transmembrane transport (II) (M191)`
## [1] 18
## 
## $`TBA (M192)`
## [1] 12
## 
## $`TBA (M193)`
## [1] 11
## 
## $`TBA (M194)`
## [1] 14
## 
## $`muscle contraction, SRF targets (M195)`
## [1] 12
## 
## $`platelet activation - actin binding (M196)`
## [1] 17
## 
## $`TBA (M197)`
## [1] 14
## 
## $`TBA (M198)`
## [1] 12
## 
## $`platelet activation & blood coagulation (M199)`
## [1] 13
## 
## $`antigen processing and presentation (M200)`
## [1] 11
## 
## $`TBA (M201)`
## [1] 18
## 
## $`enriched in extracellular matrix & associated proteins (M202)`
## [1] 12
## 
## $`TBA (M203)`
## [1] 13
## 
## $`chaperonin mediated protein folding (I) (M204.0)`
## [1] 14
## 
## $`chaperonin mediated protein folding (II) (M204.1)`
## [1] 10
## 
## $`TBA (M205)`
## [1] 11
## 
## $`Wnt signaling pathway (M206)`
## [1] 14
## 
## $`TBA (M207)`
## [1] 10
## 
## $`TBA (M208)`
## [1] 10
## 
## $`lysosome (M209)`
## [1] 10
## 
## $`extracellular matrix, collagen (M210)`
## [1] 32
## 
## $`TBA (M211)`
## [1] 11
## 
## $`purine nucleotide biosynthesis (M212)`
## [1] 12
## 
## $`regulation of transcription, transcription factors (M213)`
## [1] 21
## 
## $`TBA (M214)`
## [1] 12
## 
## $`small GTPase mediated signal transduction (M215)`
## [1] 16
## 
## $`respiratory electron transport chain (mitochondrion) (M216)`
## [1] 12
## 
## $`TBA (source: B cells) (M217)`
## [1] 10
## 
## $`TBA (M218)`
## [1] 15
## 
## $`respiratory electron transport chain (mitochondrion) (M219)`
## [1] 18
## 
## $`TBA (M220)`
## [1] 10
## 
## $`TBA (M221)`
## [1] 16
## 
## $`heme biosynthesis (II) (M222)`
## [1] 13
## 
## $`enriched in T cells (II) (M223)`
## [1] 13
## 
## $`transmembrane and ion transporters (II) (M224)`
## [1] 11
## 
## $`metabolism of steroids (M225)`
## [1] 10
## 
## $`proteasome (M226)`
## [1] 12
## 
## $`translation initiation (M227)`
## [1] 10
## 
## $`olfactory receptors (M228)`
## [1] 13
## 
## $`TBA (M229)`
## [1] 11
## 
## $`cell cycle, mitotic phase (M230)`
## [1] 12
## 
## $`respiratory electron transport chain (mitochondrion) (M231)`
## [1] 11
## 
## $`enriched for TF motif TNCATNTCCYR (M232)`
## [1] 13
## 
## $`TBA (M233)`
## [1] 15
## 
## $`transcription elongation, RNA polymerase II (M234)`
## [1] 14
## 
## $`mitochondrial cluster (M235)`
## [1] 19
## 
## $`TBA (M236)`
## [1] 14
## 
## $`golgi membrane (II) (M237)`
## [1] 17
## 
## $`respiratory electron transport chain (mitochondrion) (M238)`
## [1] 17
## 
## $`enriched in calcium signaling proteins (M239)`
## [1] 15
## 
## $`chromosome Y linked (M240)`
## [1] 17
## 
## $`TBA (M241)`
## [1] 10
## 
## $`TBA (M242)`
## [1] 11
## 
## $`TBA (M243)`
## [1] 12
## 
## $`TBA (M244)`
## [1] 13
## 
## $`translation initiation factor 3 complex (M245)`
## [1] 13
## 
## $`TBA (M246)`
## [1] 16
## 
## $`enriched in nuclear pore complex interacting proteins (M247)`
## [1] 14
## 
## $`TBA (M248)`
## [1] 11
## 
## $`TBA (M249)`
## [1] 10
## 
## $`spliceosome (M250)`
## [1] 12
## 
## $`T cell surface signature (S0)`
## [1] 27
## 
## $`NK cell surface signature (S1)`
## [1] 48
## 
## $`B cell surface signature (S2)`
## [1] 169
## 
## $`Plasma cell surface signature (S3)`
## [1] 24
## 
## $`Monocyte surface signature (S4)`
## [1] 94
## 
## $`DC surface signature (S5)`
## [1] 82
## 
## $`CD4 T cell surface signature Th1-stimulated (S6)`
## [1] 9
## 
## $`CD4 T cell surface signature Th2-stimulated (S7)`
## [1] 13
## 
## $`Naive B cell surface signature (S8)`
## [1] 54
## 
## $`Memory B cell surface signature (S9)`
## [1] 37
## 
## $`Resting dendritic cell surface signature (S10)`
## [1] 75
## 
## $`Activated (LPS) dendritic cell surface signature (S11)`
## [1] 37
# Get all unique genes
all_genes <- sort(unique(unlist(btm)))

# Create the binary matrix
binary_matrix_btm <- sapply(all_genes, function(g) {
    as.numeric(sapply(btm, function(x) g %in% x))
})

# Transpose and convert to a data frame for better visualization
#binary_matrix_btm <- as.data.frame(t(binary_matrix_btm))
binary_matrix_btm[1:5, 1:5]
##      A1CF A2BP1 A2M ABCA1 ABCA13
## [1,]    0     0   0     0      0
## [2,]    0     0   0     0      0
## [3,]    0     0   0     0      0
## [4,]    0     0   0     0      0
## [5,]    0     0   0     0      0
# Add row names
rownames(binary_matrix_btm) <- names(btm)

# Display the binary matrix
binary_matrix_btm[1:5, 1:10]
##                                                A1CF A2BP1 A2M ABCA1 ABCA13
## targets of FOSL1/2 (M0)                           0     0   0     0      0
## integrin cell surface interactions (I) (M1.0)     0     0   0     0      0
## integrin cell surface interactions (II) (M1.1)    0     0   0     0      0
## extracellular matrix (I) (M2.0)                   0     0   0     0      0
## extracellular matrix (II) (M2.1)                  0     0   0     0      0
##                                                ABCA6 ABCA8 ABCB4 ABCB6 ABCC13
## targets of FOSL1/2 (M0)                            0     0     0     0      0
## integrin cell surface interactions (I) (M1.0)      0     0     0     0      0
## integrin cell surface interactions (II) (M1.1)     0     0     0     0      0
## extracellular matrix (I) (M2.0)                    0     0     0     0      0
## extracellular matrix (II) (M2.1)                   0     0     0     0      0
dim(binary_matrix)
## [1]   50 4384
enrichment4.parametric.pos <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = 15,
                                             feature.sets = as.matrix(binary_matrix_btm),
                                             sign = "positive",
                                             statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.
## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 124 
## Set statistic: mean.diff 
## Statistical test: parametric
## Subsetting weights with positive sign
## 
plot_enrichment(enrichment4.parametric.pos, 
                factor = 'Factor15', 
                max.pathways = 15
)

plot_enrichment_detailed(enrichment4.parametric.pos, 
                         factor = 'Factor15', 
                         max.genes = 8, 
                         max.pathways = 5
)

enrichment4.parametric.neg <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = 15,
                                             feature.sets = as.matrix(binary_matrix_btm),
                                             sign = "negative",
                                             statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.
## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 124 
## Set statistic: mean.diff 
## Statistical test: parametric
## Subsetting weights with negative sign
## 
plot_enrichment(enrichment4.parametric.neg, 
                factor = 'Factor15', 
                max.pathways = 15
)

plot_enrichment_detailed(enrichment4.parametric.neg, 
                         factor = 'Factor15', 
                         max.genes = 8, 
                         max.pathways = 5
)
## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

other modules

enrichment.parametric <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = c(1:15),
                                             feature.sets = as.matrix(binary_matrix_btm),
                                             sign = "positive",
                                             statistical.test = "parametric"
)
## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.
## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 124 
## Set statistic: mean.diff 
## Statistical test: parametric
## Subsetting weights with positive sign
## 
weights <- get_weights(MOFAobject3, views = "all", factors = "Factor15", scale = F)
weights.scale  <- get_weights(MOFAobject3, views = "geneExp", factors = "Factor15", scale = T)

a= weights$ab
ggplot(a, aes(x=Factor15)) +
    geom_density() 

a['IgG_PT', ]
## [1] -0.001356319
#weights.scale$geneExp[c('CCL3', 'CXCL8', 'IL1B', 'CCL4'), ]
a = a[order(a[,1]), ]
head(a)
##     IgG1_FHA      IgG1_TT      IgG1_DT      IgG4_DT     IgG2_OVA  IgG1_FIM2/3 
## -0.009453229 -0.008087487 -0.006169425 -0.005323671 -0.005263566 -0.005103896
tail(a)
##       IgG3_TT   IgG3_FIM2/3   IgG2_FIM2/3       IgG2_TT      IgG3_OVA 
## -1.423650e-05  5.785488e-06  2.256356e-05  1.452136e-04  6.736243e-04 
##       IgG3_PT 
##  2.081636e-03
d = as.data.frame(tail(names(a), n=200))
#write_tsv(d, file='data/CCL3_MOFA_F3.tsv', col_names = T)
d
##    tail(names(a), n = 200)
## 1                 IgG1_FHA
## 2                  IgG1_TT
## 3                  IgG1_DT
## 4                  IgG4_DT
## 5                 IgG2_OVA
## 6              IgG1_FIM2/3
## 7                 IgG4_FHA
## 8                  IgG2_DT
## 9                 IgG1_PRN
## 10                 IgG_PRN
## 11                 IgG2_PT
## 12                 IgG1_PT
## 13                IgG4_OVA
## 14                IgG3_FHA
## 15                IgG2_FHA
## 16                  IgG_PT
## 17                 IgG4_TT
## 18                 IgG3_DT
## 19                IgG1_OVA
## 20                IgG3_PRN
## 21                IgG4_PRN
## 22             IgG4_FIM2/3
## 23                 IgG_FHA
## 24                 IgG4_PT
## 25                IgG2_PRN
## 26                 IgG3_TT
## 27             IgG3_FIM2/3
## 28             IgG2_FIM2/3
## 29                 IgG2_TT
## 30                IgG3_OVA
## 31                 IgG3_PT

cytokineLs

plot_top_weights(MOFAobject3,
                 view = "cytokineL",
                 factor = 15,
                 nfeatures = 5
)

plot_top_weights(MOFAobject3,
                 view = "cell_freq",
                 factor = 15,
                 nfeatures = 5
)

plot_data_scatter(MOFAobject,
                  view = "ab",         # view of interest
                  factor = 15,             # factor of interest
                  features = 5,           # number of features to plot (they are selected by weight)
                  add_lm = TRUE,         # add linear regression
                  color_by = "Meta.infancy_vac"
)

end

Sys.Date()
## [1] "2024-12-13"
getwd()
## [1] "/Users/nthrupp/Documents/Projects/CMI-PB3/models/MOFA_6assays_TP0and1_factorsOnly_split2/MOFA"
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] MOFAdata_1.18.0  MOFA2_1.12.1     lubridate_1.9.3  forcats_1.0.0   
##  [5] stringr_1.5.1    dplyr_1.1.4      purrr_1.0.2      readr_2.1.5     
##  [9] tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1    tidyverse_2.0.0 
## [13] BiocStyle_2.30.0
## 
## loaded via a namespace (and not attached):
##   [1] rlang_1.1.4           magrittr_2.0.3        matrixStats_1.4.1    
##   [4] compiler_4.3.1        mgcv_1.9-1            dir.expiry_1.10.0    
##   [7] png_0.1-8             vctrs_0.6.5           reshape2_1.4.4       
##  [10] pkgconfig_2.0.3       crayon_1.5.3          fastmap_1.2.0        
##  [13] backports_1.5.0       magick_2.8.4          XVector_0.42.0       
##  [16] labeling_0.4.3        utf8_1.2.4            rmarkdown_2.27       
##  [19] tzdb_0.4.0            bit_4.0.5             tinytex_0.52         
##  [22] xfun_0.47             zlibbioc_1.48.2       cachem_1.1.0         
##  [25] jsonlite_1.8.8        highr_0.11            rhdf5filters_1.14.1  
##  [28] DelayedArray_0.28.0   BiocParallel_1.36.0   Rhdf5lib_1.24.2      
##  [31] broom_1.0.6           parallel_4.3.1        R6_2.5.1             
##  [34] bslib_0.8.0           stringi_1.8.4         RColorBrewer_1.1-3   
##  [37] reticulate_1.39.0     car_3.1-2             jquerylib_0.1.4      
##  [40] Rcpp_1.0.13           bookdown_0.40         knitr_1.48           
##  [43] IRanges_2.36.0        Matrix_1.6-5          splines_4.3.1        
##  [46] timechange_0.3.0      tidyselect_1.2.1      rstudioapi_0.16.0    
##  [49] abind_1.4-8           yaml_2.3.10           codetools_0.2-20     
##  [52] lattice_0.22-6        plyr_1.8.9            basilisk.utils_1.14.1
##  [55] withr_3.0.1           evaluate_1.0.0        Rtsne_0.17           
##  [58] pillar_1.9.0          BiocManager_1.30.23   ggpubr_0.6.0         
##  [61] filelock_1.0.3        MatrixGenerics_1.14.0 carData_3.0-5        
##  [64] corrplot_0.92         stats4_4.3.1          generics_0.1.3       
##  [67] vroom_1.6.5           S4Vectors_0.40.2      hms_1.1.3            
##  [70] munsell_0.5.1         scales_1.3.0          glue_1.7.0           
##  [73] pheatmap_1.0.12       tools_4.3.1           data.table_1.16.0    
##  [76] fgsea_1.28.0          ggsignif_0.6.4        fastmatch_1.1-4      
##  [79] cowplot_1.1.3         rhdf5_2.46.1          grid_4.3.1           
##  [82] colorspace_2.1-1      nlme_3.1-166          basilisk_1.14.3      
##  [85] HDF5Array_1.30.1      cli_3.6.3             fansi_1.0.6          
##  [88] S4Arrays_1.2.1        uwot_0.2.2            gtable_0.3.5         
##  [91] rstatix_0.7.2         sass_0.4.9            digest_0.6.36        
##  [94] BiocGenerics_0.48.1   SparseArray_1.2.4     ggrepel_0.9.6        
##  [97] farver_2.1.2          htmltools_0.5.8.1     lifecycle_1.0.4      
## [100] bit64_4.0.5